#include <cblas.h>
#include "MatrixMultiply.h"
#include "Matrix.h"

/*
MatrixMultiply::MatrixMultiply(void) : 	m_matA(1, 1), m_matB(1, 1), m_matC(1, 1), n(1), m_lda(1), m_ldb(1), m_ldc(1)
{}
*/

MatrixMultiply::MatrixMultiply(Matrix *matA, Matrix *matB, Matrix *matC) : m_matA(*matA), m_matB(*matB), m_matC(*matC), n(matA->GetRowCount()), m_lda(matA->GetLD()), m_ldb(matB->GetLD()), m_ldc(matC->GetLD())
{
	m_pMatA = matA->GetMatrix();
	m_pMatB = matB->GetMatrix();
	m_pMatC = matC->GetMatrix();
}

void MatrixMultiply::calcNaive()
{
	for(int x = 0; x < n; x++)
	{
		for(int y = 0; y < n; y++)
		{
			double s = 0.0;
			for(int k = 0; k < n; k++)
				s = s + m_pMatA[y+(k*m_lda)] * m_pMatB[k+(x*m_ldb)];

			m_pMatC[y + x*m_ldc] = s;
		}	
	}

}
/*
void MatrixMultiply::calcNaive(int n, double *matA, double *matB, double *matC, int lda, int ldb, int ldc)
{
	for(int x = 0; x < n; x++)
	{
		for(int y = 0; y < n; y++)
		{
			double s = 0.0;
			for(int k = 0; k < n; k++)
				s = s + matA[y+(k*lda)] * matB[k+(x*ldb)];
				s = s + (*matA + k + x * lda) * (*matB + k * ldb + y);
			*matC = s; 

			matC[y + x*ldc] = s;
		}	
	}
}
*/
void MatrixMultiply::cblas_myDdot()
{
	const int eins = 1;
	const int n = m_matA.GetRowCount();
	for(int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
			m_pMatC[i+j*m_ldc] = cblas_ddot(n,&m_pMatA[i],m_lda,&m_pMatB[j*m_ldb],eins);
	}
}

void MatrixMultiply::cblas_myDaxpy()
{
	const int eins = 1;
	const int n = m_matA.GetRowCount();
	for(int j = 0; j < n; j++)
	{
		for(int i =  0; i < n; i++)
		{
			m_pMatC[i+j*m_ldc]=0.0;
		}
		for(int k = 0; k < n; k++)
		{
			double alpha = m_pMatB[k+j*m_ldb];
			cblas_daxpy(n, alpha, &m_pMatA[k*m_lda], eins, &m_pMatC[j*m_ldc], eins);
		}
	}
}

void MatrixMultiply::cblas_myDGEMV()
{
	const int eins = 1;
	const int n = m_matA.GetRowCount();
	double alpha = 1.0;
	double beta = 0.0;
	for(int j = 0; j < n; ++j)
		cblas_dgemv(CblasColMajor, CblasNoTrans, n, n, alpha, &m_pMatA[0], m_lda, &m_pMatB[j*m_ldb], eins, beta, &m_pMatC[j*m_ldc], eins);
}

void MatrixMultiply::cblas_myDGEMM()
{
	double alpha = 1.0;
	double beta = 0.0;
	const int n = m_matA.GetRowCount();
	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha, &m_pMatA[0], m_lda, &m_pMatB[0], m_ldb, beta, &m_pMatC[0], m_ldc);
}
MatrixMultiply::~MatrixMultiply(void)
{
}
